#ifndef _NMOVE
#define _NMOVE 1
#endif
template<class Density, class Proposal>//Proposal avec r et d
class Kernel_MHID{
	public:
		Kernel_MHID(Density *d, Proposal *p)
		{
			_d=d;
			_p=p;
			boost::random::uniform_01<> U;
			runif=new RandomG::Random<boost::random::uniform_01<> >(U);
			_Sig=0.8;
			_nb=0;
		}
		void Move(mat *X,int _n)//Move matrix convert to multidimentional theta 
		{		//, moves each column dimentionis in line
			int n=X->n_rows;//Y and X as pointers
			int m=X->n_cols;
			mat Y(n,m);
			Y.fill(0);
			mat xt;

			_p->Set_n(m);
			_sum=0;
			_i=0;
			for(int i=0;i<n;i++)
			{
				Y(i,span::all)=(_p->r(_i)).t();//Y[i]=p->r(X[i])
				mat xt=Y(i,span::all);
				xt=Accept(X->row(i),Y.row(i),_n);
				(*X)(i,span::all)=xt;	
				_i++;			
			}
			_accept=(double)_sum/n;

			cout << "accept = " << _sum/n;
			//////////
			//delete[] y;
		}
		void testMove(mat *X,int _n)//Move matrix convert to multidimentional theta 
		{		//, moves each column dimentionis in line
			int n=X->n_rows;//Y and X as pointers
			int m=X->n_cols;
			if(_bool==1){
				Moved.zeros(); // idicates if ith particle has moved
			}
			mat Y(n,m);
			Y.fill(0);
			mat xt;
			_p->Set_n(m);
			_sum=0;
			_i=0;
			for(int i=0;i<n;i++)
			{
				Y(i,span::all)=(_p->r(i)).t();//Y[i]=p->r(X[i])
				mat xt=Y(i,span::all);
				xt=Accept(X->row(i),Y.row(i),_n);
				(*X)(i,span::all)=xt;			
				_i++;	
			}
			_accept=(double)_sum/n;

			cout << "accept = " << _sum/n;
			//////////
			//delete[] y;
		}
		////////////////////////////////////////////////////////////////
		void Move2(mat *X, int _n)
		{
			int n=X->n_rows;//Y and X as pointers
			int m=X->n_cols;
			mat Y(n,m);
			Y.fill(0);
			mat xt;
			_p->Set_n(m);
			_sum=0;
			_i=0;
			for(int i=0;i<n;i++)
			{
				Y(i,span::all)=Q_gamma((*X)(i,span::all).t()).t();//Y[i]=p->r(X[i])
				mat xt=Y(i,span::all);
				xt=Accept2(X->row(i),Y.row(i),_n);
				(*X)(i,span::all)=xt;				
				_i++;
			}
			_accept=(double)_sum/n;

			cout << "accept = " << _sum/n;

		}
		mat Accept2(mat X, mat Y,int n){
	
			mat res=X;
			if(_bool==1)
			{
				double rho=0;
				double temp=_d->Density_utn(Y.t(),n);
				double temp2=(_d->Density_utn(X.t(),n));
				rho=(temp)-temp2;
				double u=log((*runif)());
				if(u<rho)
				{
					Dens(_i,0)=temp;
					res=Y;
					_sum+=1;
					Moved(_i,0)+=1;
				}else{
					Dens(_i,0)=temp2;
				}
			}else{
				
				double rho=0;
				double temp=_d->Density_utn(Y.t(),n);
				rho=temp-(Dens(_i,0));
				double u=log((*runif)());
				if(u<rho)
				{
					res=Y;
					_sum+=1;
					Dens(_i,0)=temp;
					Moved(_i,0)+=1;
				}

			}
			return res;
		}
		mat Q_gamma(mat gammat)
		{
			int p=gammat.n_elem;
			int d=sum(sum(gammat));
			int dc=p-d;
			 
			double un=(*runif)(); 
		 	if((dc==0) | (dc==p) | un<1)
			{	
				Distribution::Uniform U1(1,0,p-1);
				int u=as_scalar(U1.r(1));
				gammat(u,0)=1-gammat(u,0);//as_scalar((*U).r(1));	
				return gammat;
			}else{
				
				Distribution::Uniform U1(1,0,d-1);
				Distribution::Uniform U2(1,0,dc-1);
				int u1=as_scalar(U1.r(1));
				int u2=as_scalar(U2.r(1));
				int b=1;
				int i=0;
				int k=-1;
				int j=-1;
				int b1=0;
				int b2=0;
				while(b)
				{
					if(gammat(i,0)==1)
					{
						k++;
					}else{
						j++;
					}
					if(k==u1)
					{
						gammat(i,0)=0;
						b1=1;
					}
					if(j==u2)
					{
						gammat(i,0)=1;
						b2=1;
					}
					i++;
					b=1-b1*b2;
					
				//	cout << gammat;
				}
			//	cout << "\\\\\\\\" << gammat;
				return gammat;
			}
		}


	
	////////////////////////////////////////////////////////////////////////	
	mat Accept(mat X, mat Y,int n){
			mat res=X;
			//cout << X.n_rows<< ";";
			//cout << Y.n_rows;
			//	cout << Y.t();
			if(_bool){
				double rho=0;
				double rho1=_d->Density_utn(Y.t(),n);
				double rho2=_d->Density_utn(X.t(),n);
				rho=rho1-log(_p->d(Y.t(),0,_i))+log(_p->d(X.t(),0,_i))-rho2;	
				
				//rho=(_d->Density_utn(Y,n)/**_p->d(X)*/)-(/*_p->d(Y)**/_d->Density_utn(X,n));
				double u=log((*runif)());
				if(u<rho)
				{
					Dens(_i,0)=rho1;
					res=Y;
					_sum+=1;
					Moved(_i,0)+=1;
				}else{
					Dens(_i,0)=rho2;
				}
			}else{
				double rho=0;
				double rho1=_d->Density_utn(Y.t(),n);
				double rho2=Dens(_i,0);
				rho=rho1-log(_p->d(Y.t(),0,_i))+log(_p->d(X.t(),0,_i))-rho2;
				//rho=(_d->Density_utn(Y,n)/**_p->d(X)*/)-(/*_p->d(Y)**/_d->Density_utn(X,n));
			//	cout << rho << "\n";
				double u=log((*runif)());
				if(u<rho)
				{
					Dens(_i,0)=rho1;
					res=Y;
					_sum+=1;
					Moved(_i,0)+=1;
				}
			}
			return res;
		}
		void Set_s(mat X){	
			mat s=_p->Var(X);
		//	cout << "s" << s;
			_p->Set_s(s);	
		}
		double Get_accept(void)
		{
			return _accept;
		}
	
	 	void Set_Sigma(double x)
		{
			_Sig=_Sig*x;
		}
		void Set_bool(int b){_bool=b;}
		void Set_n(int b){
			_n=b;
			mat De(_n,1);
			mat M(_n,1);
			Moved=M;
			Dens=De;
		}
		double nMove(void)
		{
			double b=1;
			for(int i=0;i<_n;i++)
			{
				if(Moved(i,0)<_NMOVE)
				{
					b=0;
				}
			}
			return b;
		}
		double nMove2(mat X)
		{
			int n=nb_unique(X);		
			if((n-_nb)<10)
			{
				_nb=n;
				return 1;
			}else{
				_nb=n;
				return 0;				
			}	
		}
	private:
		Density *_d;
		mat Moved;
		mat Dens; 
		Proposal *_p;
		RandomG::Random<boost::random::uniform_01<> > *runif;
		int _n;
		int _i;
		int _nb;
		int _bool;
		double _sum;
		double _accept;
		double _Sig;
};

